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ABSTRACT 

Is the turbulence in cluster-forming regions internally driven by stellar outflows or the consequence of 
a large-scale turbulent cascade? We address this question by studying the turbulent energy spectrum 
in NGC 1333. Using synthetic ^'^CO maps computed with a snapshot of a supersonic turbulence 
simulation, we show that the VCS method of Lazarian and Pogosyan provides an accurate estimate 
of the turbulent energy spectrum. We then apply this method to the ^"^CO map of NGC 1333 from 
the COMPLETE database. We find the turbulent energy spectrum is a power law, E{k) oc k~f^ , in 
the range of scales 0.06 pc < € < 1.5 pc, with slope (3 — 1.85 ± 0.04. The estimated energy injection 
scale of stellar outflows in NGC 1333 is £inj ~ 0.3 pc, well resolved by the observations. There is 
no evidence of the flattening of the energy spectrum above the scale ^inj predicted by outflow-driven 
simulations and analytical models. The power spectrum of integrated intensity is also a nearly perfect 
power law in the range of scales 0.16 pc< I <7.9 pc, with no feature above ^in^y We conclude that the 
observed turbulence in NGC 1333 does not appear to be driven primarily by stellar outflows. 
Subject headings: ISM: kinematics and dynamics — stars: formation — turbulence 



1. INTRODUCTION 

The structure of density and velocity fields in gi- 
ant molecular clouds can be characterized by extended 
power laws (e.g. Heyer and Brunt 2004; Padoan et al. 
2004, 2006), with scaling exponents consistent with those 
of numerical simulations of supersonic turbulence (e.g. 
Padoan et al. 2007; Kritsuk et al. 2007, 2009a,b). The 
existence of these scaling laws and of one-point statistics 
of turbulent flows (for example the probability distribu- 
tion of gas density) is an important assumption in statis- 
tical theories of star formation aimed at the prediction of 
the stellar initial mass function (Padoan and Nordlund 
2002, 2004; Hannebelle and Chabrier 2008, 2009) and the 
star formation rate (Krumholz and McKee 2005; Padoan 
and Nordlund 2009). 

Is this assumption of a universal large-scale turbulent 
cascade valid in cluster- forming regions, believed to be 
the birthplace of most stars, or is the turbulence there 
driven internally by stellar outflows? As an alternative 
to the idea that stellar clusters are formed in one cross- 
ing time (Elmegreen 2000), Tan, Krumholz, and Mc- 
Kee (2006) propose that stars are formed in protocluster 
clumps over several dynamical times, while the turbu- 
lence is driven by stellar outflows. This scenario is inves- 
tigated analytically by Matzner (2007) and numerically 
by Nakamura and Li (2007), CarroU et al. (2007), and 
Wang et al. (2009). 

Nakamura and Li (2007) find that the velocity power 
spectrum in their outflow-driven simulation is very shal- 
low above a characteristic scale of energy injection by 
outflows, liny They derive ^inj ~ 0.3 pc, assuming their 
computational domain has a size L = 1.5 pc and a total 
mass TOtot = 929 mQ , characteristic of cluster- forming re- 
gions. This value of ^inj and the flattening of the power 
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spectrum at larger scales is consistent with the analyti- 
cal resuhs in Matzner (2007). Carroll et al. (2009) con- 
firm the results of Nakamura and Li (2007) at higher nu- 
merical resolution, producing velocity power spectra even 
more clearly peaked at the outflow driving scale. With 
the values of mass, column density, and outflow momen- 
tum per unit mass, Vc ~ 50 km/s, adopted by Naka- 
mura and Li (2007), equation (38) in Matzner (2007) 
gives €inj = 0.27 pc. The injection scale is only weakly 
dependent on v^, mtot, and the star formation rate, and 
is approximately given by ^inj « 0.079 Scgs^^^* pc, where 
Scgs is the column density in g/cm^. Our conclusions do 
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not critically depend on a precise knowledge of 
long as scales just above £inj are resolved in the observa- 
tions. 

The prototype of outflow-driven cluster-forming re- 
gions chosen by Nakamura and Li (2007), Matzner 
(2007), and Carroll et al. (2009) is NGC 1333 in Perseus. 
In this Letter, we derive the power spectrum of integrated 
intensity, P/(fc), and velocity, Pv{k), in NGC 1333, based 
on the data from the Five College Radio Astronomy 
Observatory (FCRAO) survey of the Perseus molecular 
cloud complex (Ridge et al. 2006), publicly available 
from the COMPLETE website. We show that the scale 
^inj is resolved by the observations, and both power spec- 
tra, Pi{k) and Py{k), are consistent with power laws ex- 
tending to large scale, with the same slope found in su- 
personic turbulence simulations and in the Perseus com- 
plex on larger scale (Padoan et al. 2006). The constant 
slopes of Pv{k) and Pi{k) above the scale £inj suggest 
that outflows may not be the dominant driving mecha- 
nism. In a recent paper appeared after the submission 
of this Letter, Brunt, Heyer, and Mac Low (2009) have 
reached the same conclusion for NGC 1333 and for other 
molecular cloud regions using a different method. 

2. POWER SPECTRUM FROM THE VCS METHOD 

We first confirm the validity of the velocity coordi- 
nate spectrum (VCS) method of Lazarian and Pogosyan 
(2006), using synthetic maps of the J=l-0 fine of ^'^CO, 
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Fig. 1. — Power spectra of integrated intensity, Pj{k), from the 
three simulated ^''CO maps. The wavenumber corresponding to 
the scale of 1 pc is marked by an arrow. The plots are compensated 
by a factor of fc, so the absolute values of the slopes given in the 
figure are the values of /By — I for the three maps. 



Fig. 2. — Velocity coordinate power spectra from the three sim- 
ulated '^^CO maps. The wavenumber above which the spectra are 
expected to steepen (in the absence of thermal broadening), fc„ 
is marked by an arrow. The absolute values of the least squares 
fit exponents given in the figure are the values of /3„ for the three 
maps. 



computed with a non-LTE radiative transfer code (Ju- 
vela and Padoan 2005), based on the density and veloc- 
ity fields from a snapshot of a simulation of supersonic 
hydrodynamic turbulence with rms Mach number M—6. 
These synthetic maps are the same used by Padoan et al. 
(2006) to test the VGA method of Lazarian and Pogosyan 
(2000) (Padoan et al. (2006) applied the VGA method 
to the full FCRAO Perseus map and found a power law 
turbulent energy spectrum, E{k) oc k'^Py{k) oc k~^, with 
the exponent /3 = 1.81 ± 0.10). 

The simulation is carried out with the Enzo code, 
developed at the Laboratory for Computational Astro- 
physics by Bryan, Norman and collaborators (Norman 
and Bryan 1999). Enzo is a public domain Eulerian grid- 
based code (see http://lca.ucsd.edu/projects/enzo) that 
adopts the Piecewise Parabolic Method (PPM) of Golella 
and Woodward (1984). We use an isothermal equation 
of state, periodic boundary conditions, initially uniform 
density and random large-scale velocity. The turbulence 
is forced in Fourier space only in the wavenumber range 
1 < fc < 2, where k — 1 corresponds to the size of 
the computational domain that contains 1, 024^ compu- 
tational zones (for details see Kritsuk et al. 2007). 

The radiative transfer calculations assume a box size 
of 5 pc, a mean density of 10^ cm^"^, a mean kinetic 
temperature of 10 K, an rms Mach number A4=6 (con- 
sistent with the turbulence simulation) and a uniform 
^■^GO abundance of 10^^. These values were chosen as 
a generic reference model, not tailored to the Perseus 
molecular cloud complex. The density and velocity data 
cubes are resampled from 1,024'^ to 256'^ zones. The re- 
sampling of the data cubes has several advantages: i) It 
yields density and velocity fields with power spectra that 
are power laws almost up to the new Nyquist frequency; 
ii) It speeds up the radiative transfer calculations; iii) 
It generates a map of synthetic spectra with a range of 
scales comparable to that of the map of NGG 1333. 

As a result of the radiative transfer calculations, we 



compute three maps (one in each of the three axis di- 
rections) of spectral profiles of the ^'^GO (J=l-0) line 
intensity, T(x, w), where x is the map position and v 
the velocity channel. We compute r(x, v) for 280 ve- 
locity channels and 256^ map positions, so we obtain a 
position-position- velocity data cube of 256 x 256 x 280 el- 
ements, for each of the three coordinate directions. The 
width of the velocity channels is Av — 0.025 km/s, giv- 
ing a total velocity range of 7 km/s across the 280 chan- 
nels. The integrated intensity, /(x), is a two-dimensional 
map obtained as the sum of the line intensity over all 
the velocity channels, multiplied by the channel width, 
/(x) = S„r(x,w) Av. 

We call Piik) the two-dimensional power spectrum of 
the integrated intensity map, /(x), and approximate it 
with a power law fit such that, Pi{k) oc k^^' . We then 
call E{k) the energy spectrum, given by E{k) oc k'^Py(k), 
where Pvik) is the three dimensional power spectrum 
of the velocity field, also assumed to be a power law, 
E{k) oc k^Pv{k) oc fc-^. The VGS method of Lazar- 
ian and Pogosyan (2006) allows the determination of 
the energy spectrum (the exponent /5), directly from 
the one dimensional power spectrum in velocity space, 
Pi(x, ky) = \ j dv e*''""^ r(x, u)!^, that is the power spec- 
trum of individual spectral line profiles. To reduce the 
effect of instrumental noise and of intrinsic variations 
in the power spectrum across the map, a global power 
spectrum, P„(fc„), is obtained as the sum of Pi(x, fci,) 
over a region of the map, Pv{kv) = SxPi(x, ky). If (3y is 
the exponent of the velocity coordinate power spectrum, 
Py{ky) OC fc"^", the turbulent energy spectrum exponent 
is given by /3 = 1 + 2/f3v, if two conditions are satisfied: 
1) a steep density power spectrum, Pp{k) oc k~^'', with 
l3p > 3, and 2) small velocity wavenumbers, k^ < kv^Ax, 
with ky^Ax = (Ax/S')-™/^, where A X is the observa- 
tional beam size, S is the size of the mapped region, and 
m is the exponent of the second order velocity structure 
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function, m — f3 ~ 1. 

The first condition, /3p > 3, is fiere approximated as 
Pi > 3, because the power spectrum of the integrated 
intensity, Pi(k), is expected to provide a good estimate 
of the density power spectrum. This is confirmed by our 
study of the synthetic spectral maps (see below) , show- 
ing that fi — \ + 2/ fiy while /?/ > 3. The second con- 
dition means that the exponent must be evaluated 
only for large velocity scales (small velocity wavenum- 
bers) relative to the characteristic turbulent velocity at 
the scale of the beam size, Aa:. We have made the ve- 
locity (hence k^) nondimensional, v v/{Ntot x Aw), 
where iVtot is the total number of velocity channels used 
to compute the power spectrum, and have defined the 
velocity wavenumber as fc„ = 1/v, so fe^, = 1 corresponds 
to the total velocity range used to compute the power 
spectrum. Because the characteristic turbulent velocity 
scales as v{i) oc i"^/"^, the velocity wavenumber scales as 
ky{£) oc explaining the dependence of ky^Ax on 

the beam size, Ax. 

Before computing the power spectra, we apply a Gaus- 
sian beam with FWllM=dx, where dx is the numerical 
mesh size, and add Gaussian noise to the synthetic spec- 
tra, to a level comparable to that of the FCRAO map 
of NGC 1333. The power spectra P/(fc) and Pv{ky) are 
shown for all three directions in Figures 1 and 2 respec- 
tively. Pi{k) is very well described by a power law in the 
approximate range 1 < fc < 50 (fc = 1 corresponds to 
the linear size of the computational domain). Averaging 
the results of the three maps, we obtain the exponent 
Pi = 3.10 ± 0.08, where the uncertainty is the standard 
deviation of the three measurements, a little larger than 
the standard deviation of each of the three individual 
least squares fits. In Padoan et al. (2006) we measured 
the same power spectrum, but without adding noise to 
the synthetic maps, and only in the direction of the x- 
axis. We obtained a slope of 2.99 ± 0.08, consistent with 
the slope of 3.07 found here for the same direction. This 
slope is also consistent with that of the density power 
spectrum measured directly from the three-dimensional 
snapshot, /3p = 3.0. 

The least squares fits to the Pv{kv) spectra are com- 
puted in the range 1 < A;„ < ky^^x, where kv,Ax = 
256'^''^^-'/^ — 10, to satisfy the second condition for the 
validity of the VCS method. As shown in Figure 2, the 
power spectra P«(fc„) are power laws up to fc„ k^^Ax, 
and become steeper at larger values of ky , as predicted by 
the theory. The steepening is here reduced by the effect 
of the added noise. Because the condition /?/ > 3 is sat- 
isfied, we derive the energy spectrum from the relation 
P =l + 2/l3y. We obtain /3 = 1.82 ± 0.04, where the un- 
certainty is the standard deviation of the exponents from 
the three maps. This result is consistent with the slope of 
1.8 derived directly from the three-dimensional velocity 
field of the simulation snapshot. The value of Pj is very 
close to the value of 3 that separates the regimes of steep 
and shallow density power spectra discussed in Lazarian 
and Pogosyan (2000, 2006). This is not a problem, be- 
cause in the case of shallow density power spectra, the 
relation between /3 and Py would be either the same or 
P = l + 2{Pp — 2)/Py, depending on some other condition. 
The two formulas give the same result for /3p ^ 3. 

We conclude that the VCS method provides a precise 
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Fig. 3. — Integrated intensity map of NGC 1333, from the region 
of the COMPLETE survey of Ridge et al. (2006) used in this 
work. The power spectra are also computed within the smaller 
square region marked by dotted lines, with essentially the same 
result. The two segments indicate the size of the outflow injection 
scale, and the computational box size, LlNi in the simulation 
of Nakamura and Li (2007), assuming a distance to NGC 1333 of 
280 pc. Assuming this distance, the map size is 7.9x7.9 pc, while 
the smaller map is 3.2x3.2 pc. 



estimate of the turbulent energy spectrum, by fitting the 
slope of Py{ky) for ky < ky^Ax- The ky range of the power 
law is limited. It is only a factor of 10 for the synthetic 
maps (and only a factor of 4 for the map of NGC 1333), 
because it is essentially given by the spatial dynamical 
range (256 for the synthetic maps) to the power m/2 — 
[P — l)/2, as explained above. However, the average of 
the power spectra over many positions on the map (256^ 
for both the synthetic maps and the map of NGC 1333), 
yields a Py{ky) spectrum that is extremely smooth and 
very well approximated by a power law, providing an 
accurate estimate of the slope of the turbulent energy 
spectrum. 

3. THE POWER SPECTRA OF NGC 1333 

We now apply the VCS method to the J=l-0 ^^CO sur- 
vey of the Perseus molecular cloud complex carried out 
with the FCRAO 14 m antenna by Ridge et al. (2006). 
The grid spacing of the survey is 23" , and the beam 
size 46". The velocity-channel size is 0.06 km/s. We se- 
lect a squared region of 100x100 arcmin, approximately 
centered on NGC 1333. This region contains 256 x 256 
spectra, like the synthetic maps. Figure 3 shows the in- 
tegrated intensity map of that region. It also shows the 
size of ^inj = 0.3 pc and Lln = 1-5 pc, the outflow in- 
jection scale and the size of the computational domain 
in Nakamura and Li (2007) and CarroU et al. (2009), 
assuming a distance to NGC 1333 of D = 280 pc. For 
this distance, the size of the map is S" = 7.9 pc, the 
grid spacing 0.03 pc, and the beam size Aa; — 0.06 pc. 
The injection scale is therefore well resolved. The un- 
certainty on the distance to NGC 1333 is not very large. 
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Fig. 4. — Power spectrum of the integrated intensity map of 
NGC 1333 shown in Figure|3] compensated by a factor k, kPj(k). 
The wavenumbers corresponding to the outflow injection scale, fci^j , 
and the computational box size, /clNi in th^ simulations of Naka- 
mura and Li (2007) and Carroll et al. (2009), assuming a distance 
to NGC 1333 of 280 pc, are marked by two arrows. The wavenum- 
ber A: = 1 corresponds to the map size of 7.9 pc. The power spec- 
trum is well appoximated by a power law in the range 1 < fc < 50, 
and the least squares fit is computed in the range 2 < fc < 40. 



Fig. 5. — Velocity coordinate power spectrum from the PPV data 
cube corresponding to the map of NGC 1333 shown in Figure |3] 
The velocity wavenumbers corresponding to the outflow injection 
scale, fc„_inj, and the computational box size, fc^i^LNi in th^ sim- 
ulations of Nakamura and Li (2007) and Carroll et al. (2009), 
assuming a distance to NGC 1333 of 280 pc, are marked by two ar- 
rows. The dynamical range in velocity coordinates is compressed, 
relative to the range in spatial scales, as explained in the text. The 
least squares fit, computed in the range 2 < fc„ < 8, gives the value 
of Pv = 2.3 ±0.04, used in the VCS method to derive the exponent 
of the turbulent energy spectrum, /3 = 1.85 ± 0.02. 



220 pc< D <350 pc (Borgman and Blaauw 1964; Her- 
big and Jones 1983; Cernis 1990). Even assuming the 
largest value oi D = 350 pc, £inj would still be 4 times 
larger than the beam size. 

We compute the power spectra Pi{k) and Pv{kv) 
within that region, without any correction for the ef- 
fects of beam and noise. This is justified because we do 
not compute the least squares fit at very large values of 
ky and k, and because we have confirmed the validity of 
the VCS method using synthetic maps where noise was 
added to a comparable level as in the map of NGC 1333. 
The Pi{k) power spectrum is shown in Figured The 
least squares fit, computed in the range 2 < fc < 40, 
gives fii = 3.19 ± 0.04. The values of the wavenumbers 
fcinj and /cln, corresponding to the scales ^inj = 0.3 pc 
and Lln = 1-5 pc, are marked by two arrows. The power 
spectrum is an almost perfect power law in the range of 
wavenumbers 1 < fc < 50, corresponding to the range of 
scales 0.16 pc< I <7.9 pc (fc = 1 corresponds to the size 
of the map), with no significant feature around ^inj- 

Figure [5] shows the velocity coordinate power spec- 
trum, Py{ky), and its least squares fit in the range 
2 < fc« < 8, giving = 2.34 ± 0.04. This is the range 
of fc„ values expected from the VCS theory, because here 
kv,Ax = 128*''"^)/^ = 7.7. As explained above, the dy- 
namical range in k^ is compressed with respect to that in 
fc, as shown by the positions of fc^^nj and fct,,LN, marked 
by two arrows. However, the measurement of the slope 
is still very accurate because the power spectrum is ex- 
tremely smooth, as the result of averaging over the whole 
map, and very well represented by a power law. The in- 
jection scale falls in the very middle of the power law 
range, at fct,_inj — 3.95, and there is no sign of a variation 
of the slope within a range of values of fc^ corresponding 
to Ax < ^ < Lln, or 0.06 pc < f < 1.5 pc. 

Using the relation (3 = 1 + 2/f3v, we obtain (3 = 
1.85 ± 0.02, where the standard deviation is now derived 
directly by that of the least squares fit for This un- 



certainty may seem unrealistically small. Using the syn- 
thetic maps we have derived a standard deviation of 0.04 
from the three values of /3. If we instead derive the un- 
certainty directly from the standard deviation of each fit 
of (3y , we find that the average of the three uncertainty is 
0.02, exactly like in the observational map. We therefore 
estimate that a more realistic uncertainty for the derived 
slope in NGC 1333 should be « 0.04. 

To further verify this result, we also select a smaller 
region of 40 x 40 arcmin (3.2 x 3.2 pc) around NGC 1333, 
shown as a dotted square in Figure O We obtain /3„ — 
2.29 ± 0.09, corresponding to /3 = 1.87 ± 0.03, consistent 
with the value estimated on the larger scale, though with 
a slightly larger uncertainty due to the smaller dynamical 
range and the smaller number of Pv{ky) spectra in the 
map. The good consistency between the two maps is 
not surprising, because the sum of the Pv{ky) spectra is 
dominated by the positions with the strongest integrated 
intensity that are mostly located within the inner map. 

We conclude that the VCS method yields a turbulent 
energy power spectrum with slope (3 = 1.85±0.04, in the 
range of scales 0.06 pc < £ < 1.5 pc, with no evidence of 
internal driving by outflows at the scale iinj = 0.3 pc. 

4. CONCLUSIONS 

It has been suggested that protocluster clumps actively 
form stars for several dynamical times, while being sup- 
ported by the turbulence internally generated by stellar 
outflows (e.g. Tan, Krumholz, and McKee 2006). This 
scenario is simulated by Nakamura and Li (2007) and 
Carroll et al. (2009), with physical parameters charac- 
teristic of NGC 1333. They find an energy injection scale, 
^inj ~ 0.3 pc, above which the turbulent energy spectrum 
is very shallow, up to the size of 1.5 pc of their computa- 
tional domain. This value of and the flattening of the 
power spectrum above that scale are consistent with the 
results of the analytical model of Matzner (2007). If the 
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turbulence in cluster-forming regions is primarily driven 
by outflows, the turbulent energy spectrum should flat- 
ten above the scale £i„j , up to a scale where the external 
turbulent cascade dominates again. 

We have searched for such a feature in the energy spec- 
trum of NGC 1333, the prototype of outflow-driven re- 
gions in the cited theoretical works. To this aim, we 
have applied the VCS method of Lazarian and Pogosyan 
(2006) to the ^^CO FCRAO map of NGC 1333, af- 
ter successfully testing the method on simulated data. 
We have found the energy spectrum is a power law, 
E{k) oc fc2p„(fc) oc fc-'', with 13 = 1.85 ± 0.04, in the 
range of scales 0.06 pc < ^ < 1.5 pc, with no evidence 
of internal driving by outflows at €inj = 0.3 pc. The 
power spectrum of integrated intensity is also a power 
law, Pi{k) oc k-f^', with f3i = 3.19 ± 0.04, in the range 
of scales 0.16 pc< f. <7.9 pc, with no significant features 
above the predicted injection scale. 

These power spectra are consistent with those of 



large-scale driven simulations of supersonic turbulence, 
and with those measured on larger scale in Perseus and 
other molecular cloud complexes. Although outflows 
from young stars are present in NGC 1333, the large- 
scale turbulent cascade appears to be the main energy 
source. The turbulence in NGC 1333 is either currently 
driven by significant mass inflow from larger scales, or 
in the process of being dissipated, until the time when 
winds, outflows, and ionizing radiation from stars will 
completely disperse the star-forming gas. 
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